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We propose a pseudopotential for the electron-electron Coulomb interaction to improve the ef¬ 
ficiency of many-body electronic structure calculations. The pseudopotential accurately replicates 
the scattering properties of the Coulomb interaction, and recovers the analytical solution for two 
electrons in a parabolic trap. A case study for the homogeneous electron gas using the diffusion 
Monte Carlo and configuration interaction methods recovers highly accurate values for the ground 
state energy, and the smoother potential reduces the computational cost by a factor of ~ 30. Finally, 
we demonstrate the use of the pseudopotential to study isolated lithium and beryllium atoms. 
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Electron-electron interactions drive chemical reactions, 
govern material properties, and conspire to form strongly 
correlated phases. Despite the widespread and im¬ 
portant consequences of electronic correlations, lead¬ 
ing computational techniques such as diffusion Monte 
Carlo (DMC) [I], truncated configuration interaction 
(Cl) mm, Mpller-Plesset theory [3], coupled cluster the¬ 
ory m, and FI2 methods [5]. These approaches are very 
expensive for real-life systems because the divergence in 
the electron-electron Coulomb interaction must be sam¬ 
pled finely mm- Here we propose a pseudopotential 
that accurately replicates the scattering properties of the 
Coulomb interaction, delivers the ground state energies 
within chemical accuracy of Ikcalmofo^, but does not 
diverge, which reduces the computational cost of both 
DMC and Cl by a factor of ~ 30. 

Pseudopotentials were first introduced by Hellmann [9] 
to describe the attractive electron-ion interaction in 
molecules and solids. Integrating out the core electrons 
that screen the central ion leaves a pseudopotential for 
the valence electrons. The reduction in the number of 
electrons and the greater smoothness of the electron-ion 
pseudopotential provides computational advantages that 
led to their widespread adoption in electronic structure 
calculations, including density functional theory m and 
DMC methods [In- 

First principles approaches must still account for the 
divergent repulsive electron-electron interaction that ne¬ 
cessitates fine sampling mm- The Kato cusp condi¬ 
tions [IMS] enforce a wavefunction with a kinetic en¬ 
ergy divergence that cancels the Coulomb divergence, 
leaving a remnant finite discontinuity in the local en¬ 
ergy which is evaluated with the elec¬ 

trons at point R in configuration space. There have been 
attempts to apply a local density solution to the short- 
ranged behavior Ennui- It was also proposed to intro¬ 
duce a soft-Coulomb operator either in real space ED , or 
reciprocal space ESI- Another attempt was to split the 
Coulomb interaction into a short and long-ranged com¬ 
ponent, so that they could be handled separately ESj- 
However, at present pseudopotentials are not generally 


used to smooth the electron-electron interaction. 

We develop an accurate electron-electron pseudopoten¬ 
tial for electrons scattering with any energy and angular 
momentum. We build on the formalism used to construct 
a pseudopotential for the contact interaction found in ul¬ 
tracold atomic gases ESI- This formalism is somewhat 
different from the standard pseudopotential approach de¬ 
veloped for attractive electron-ion interactions that fo¬ 
cuses on discrete bound state energies ESHIH], although 
it can be extended to scattering states EO] - The proposed 
pseudopotential is identical to the Coulomb interaction 
outside of a cut-off radius where many-body physics be¬ 
comes important. The pseudopotential delivers all of the 
physics of the Coulomb interaction but does not diverge, 
so that the ground state can be determined efficiently. 
After developing the pseudopotential in the two-body 
scattering problem, we test it on the analytically solv¬ 
able system of two electrons in a parabolic trap [3T| . 

We study the applicability, accuracy, and portability 
of the pseudopotential for a homogeneous electron gas 
(HEC) using two methods: DMC in which the use of 
the pseudopotential reduces the required time-step, and 
Cl in which the pseudopotential reduces the size of the 
plane-wave basis set required. The pseudopotential de¬ 
livers chemical accuracy, and at the same time reduces 
the computational cost of both techniques by a factor of 
~ 30. Finally, we test the pseudopotential on two inho¬ 
mogeneous systems, the isolated lithium and beryllium 
atoms. 


CONSTRUCTION OF THE PSEUDOPOTENTIAL 

To construct the pseudopotential, we adopt the for¬ 
malism of Ref. ES] and study the two-body problem: two 
electrons in their center-of-mass frame with wave vector 
k > 0 and angular momentum quantum number £. The 
Hamiltonian in atomic units is ~^V” + 
V{r)ip = and the repulsive Coulomb potential is 
V{r) = Ijr. The proposed pseudopotential is identical 
to the Coulomb potential outside of a cutoff radius c. 
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(a) Pseudopotential (b) Error with cutoff 




FIG. 1. (Color online) (a) The interaction potentials: the 
Coulomb potential is shown in red, the pseudopotential with 
cutoff radius c = ro is in cyan, and the pseudopotential with 
cutoff radius c = 2ro is in blue, (b) The error in the loga¬ 
rithmic derivative of the scattering wavefunction with cutoff 
radius for an electron gas with — 2. S shows the error 
summed over angular momentum channels, 5o is the contri¬ 
bution from the i = 0 channel and di from the i = 1 chan¬ 
nel. (c) The pseudopotential error for a range of values, 
(d) The error in the logarithmic derivative of the scattering 
wavefunction with incident wave vector for the £ = 0 and 
£ = 1 scattering channels. The filled blue curve plotted on 
an arbitrary linear scale on the secondary y-axis shows the 
weighting factor pe{k) used in evaluating the overall error in 
the logarithmic derivative. 


and at the cutoff it is continuous and has a continuous 
first derivative. At small electron-electron separation r, 
the pseudopotential can be chosen to be softer than the 
Coulomb interaction so that on electron coalescence at 
r = 0 it is finite and has zero gradient to remove pos¬ 
sible divergences and discontinuities in the local energy, 
thereby reducing the variance in our estimate of the total 
energy. These considerations suggest a pseudopotential 
of the form 





c 


r > c, 
r 

with variational freedom introduced by a polynomial ex¬ 
pansion of order = 6. To determine the parame¬ 
ters {vi} we calculate the scattering states. The scat¬ 


tering states for the Coulomb interaction can be 

solved exactly in terms of Whittaker functions, whereas 
the scattering solution (j)k,eir) from the pseudopotential 
is solved numerically. The difference between the scat¬ 
tering properties of the two potentials is characterized by 
the mean square error in the logarithmic derivative of the 
scattering wavefunction at the cutoff radius 
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which is summed over all angular momentum channels 
£ = {0, ...,6} and integrated over all possible scat¬ 
tering wave vectors 0 < A: < fcp encountered in an 
electron gas with Fermi momentum fcp |29j . Following 
Ref. [5^, we weight the importance of different scatter¬ 
ing states by a factor pe{k), which is chosen to repli¬ 
cate the density of scattering states in a Hartree-Fock 
trial wavefunction for a homogeneous electron gas where 
p^(k) = f dqnp(q)np(k -f q)/\/ (2£ -|- 1)!!, and np is the 
Fermi-Dirac distribution function. We select the varia¬ 
tional parameters {u© that minimize 6^, which gives a 
pseudopotential whose scattering closely replicates the 
Coulomb interaction. We associate the length scale 
ro = (97r/4)^/^/fcp with the typical electron separation 
and characterize an electron gas with the standard den¬ 
sity parameter r^ = ro/ap, where op is the electron Bohr 
radius) 

In Fig. [^a) we examine two of the pseudopotentials 
constructed to be used in an electron gas with density 
Ts = 2 . At small r the pseudopotential is flat to ensure 
that the wavefunction is smooth. The pseudopotential 
is therefore weaker than the Coulomb potential but, to 
give the same net scattering strength, the pseudopoten¬ 
tial must exceed the Coulomb potential at intermediate 
r, before they merge at the cutoff radius. The figure also 
shows that on reducing the cutoff radius the pseudopo¬ 
tential approaches the Coulomb potential. Therefore, the 
pseudopotential should recover the scattering properties 
of the Coulomb potential with increasing accuracy as the 
pseudopotential cutoff radius is reduced. We verify this 
in Fig. Sb) where the error in the logarithmic derivative 
of the scattering wavefunction falls with cutoff radius as 
~ (c/ro)^'®. The £ = 0 and £ = 1 channels provide 
similar contributions to the error in the scattering wave 
function. Now that we have tested the pseudopotential 
developed for an electron gas at rg = 2, we develop and 
test pseudopotentials to be applied to electron gases with 
the full range of densities 1 < Tg < 16 that can be found 
in real-life systems. Fig. [^c) shows that the average er¬ 
ror in the logarithmic derivative S ^ 10“"^ is small com¬ 
pared with the typical scattering phase shift 27r over a 
wide range of electron gas densities, demonstrating that 
the pseudopotential accurately reproduces the two-body 
scattering properties of the Coulomb interaction. 

The error in the logarithmic derivative of the wavefunc- 
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(a) ^ = 0 wavefunction (b) Two-body scattering 




FIG. 2. (Color online) (a) Wavefunction of the relative mo¬ 
tion of two electrons in a harmonic trap in the i = Q angular 
momentum state. The green curve shows the wavefunction for 
non-interacting electrons, the red shows electrons interacting 
via the Coulomb interaction, and blue shows interactions via 
the pseudopotential, (b) The error per electron in the esti¬ 
mate of the ground state energy of two opposite-spin {£ = 0) 
and same-spin (^ = 1) electrons in a parabolic trap as a func¬ 
tion of cutoff radius. The vertical blue dotted line shows the 
typical separation of the electrons in the harmonic trap. 

tion averaged over the incident wave vectors of electrons 
scattering off the pseudopotential is small. To under¬ 
stand how this is achieved, we examine in Fig. [^d) the 
variation of the error in the logarithmic derivative with 
respect to the incident wave vector. The i = 0 chan¬ 
nel has a quadratic form that crosses zero error twice, 
whereas the I = \ channel has an error that crosses zero 
only once. The variational freedom in the pseudopoten¬ 
tial has been used to minimize the error around fc « O.Sfcp 
where the density of scattering states is largest, sacrific¬ 
ing accuracy at higher incident wave vectors. 

With the pseudopotential providing phase shifts with 
an error of only ^ 10““^, we are well-positioned to test its 
performance in a many-body setting. We first study an 
idealized system with an analytical solution to provide 
an exact benchmark: two electrons in a parabolic trap. 
We also study systems that cannot be solved analyti¬ 
cally: the HEG with two complementary methods; DMC 
and Cl; and we also study isolated lithium and beryllium 
atoms. This allows us to assess the performance and ac¬ 
curacy of the pseudopotential, and verify its portability. 


TWO ELECTRONS IN A PARABOLIC TRAP 

Now that we have constructed the Coulomb pseudopo¬ 
tential and calibrated it against the phase shift of two 
atoms scattering in a vacuum, we evaluate the accuracy 
of the pseudopotential in a second analytically soluble 
system: Hooke’s atom, two interacting electrons trapped 
in the parabolic well ma;^r^/2 [3T]. This problem re¬ 
ceived early numerical attention |32H34) . and was more 
recently studied with coupled cluster methods [351 ES] ■ 


We solve separately for opposite- and same-spin electrons 
as the relative wavefunctions differ due to fermion anti¬ 
symmetry. 

We solve for the energy of two interacting electrons 
in the parabolic trap within the center-of-mass frame 
in which the interacting Hamiltonian for relative mo¬ 
tion is H A) + ccVVd + ^(£+ l)/r2 + I/(r), 

where V{r) = 1/r is the Coulomb interaction in atomic 
units. For the special case of w = 1/8 this model can 
be solved analytically for the i = 0 (opposite-spin elec¬ 
trons) ground state giving eigenenergy E = 5/4 (the non¬ 
interacting center-of-mass Hamiltonian has energy 3/4 
giving a total energy E = 2). On replacing the inter¬ 
action potential by a pseudopotential, the Hamiltonian 
for relative motion can be solved numerically and the 
ground state energy compared with the exact solution for 
the Coulomb interaction. When constructing the pseu¬ 
dopotential we chose a maximum energy of the scattering 
states that we integrate over in Eqn. ([^. We take this 
to be the energy per electron in the interacting system, 
E = 1. 

The parabolic trap is an ideal setting to compare the 
ground state wavefunction predicted by the Coulomb 
interaction with that from the pseudopotential. In 
Eig.j^a) we show the £ = 0 (i.e., opposite-spin electrons) 
ground state wavefunction for relative electron motion. 
Eirstly, to orient the discussion we show the wavefunc¬ 
tion for non-interacting electrons, which is a Gaussian 
that is smooth at electron coalescence. The wavefunction 
for the Coulomb interaction has a gradient discontinuity 
at electron coalescence which provides a divergent kinetic 
energy that cancels the divergence in the Coulomb inter¬ 
action. In general the gradient discontinuity is difficult 
to capture numerically and it hinders computational ap¬ 
proaches. However, the smooth pseudopotential provides 
a wavefunction that is smooth over all space including 
at electron coalescence, which should aid computational 
methods. 

In Fig. I^b) we study the error in the £ = 0 ground 
state energy when varying the cutoff radius, which is the 
control parameter for adjusting the accuracy of the pseu¬ 
dopotential. The error in the ground state energy with 
the cutoff set to the typical electron separation, 1/Vw, 
is 2 X 10“® au per electron. With decreasing cutoff ra¬ 
dius c the pseudopotential approaches the Coulomb in¬ 
teraction and the accuracy further increases, varying as 
~ . This scaling in error with cutoff radius is 

similar to that seen in the error in the logarithmic deriva¬ 
tive of the scattering wavefunction shown in Fig. [^b), 
which varies as ^ (c/ro)^'®. 

The interactions between opposite-spin and same-spin 
electrons both make important contributions to the total 
energy in many systems. Therefore, we next study the 
ground state energy of same-spin electrons in a parabolic 
trap. This requires a spatially anti-symmetric ground 
state, and so we require the system with £ = 1. Here 
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the system is analytically soluble with w = 1/16 giving 
an energy oi E = 21/16. In Fig. [^b) we study the er¬ 
ror in the prediction of the i = 1 energy. The error in 
the ground state energy with the cutoff set to the typical 
electron separation, l/-\/w, is 3 x 10“® au per electron. 
With decreasing cutoff radius c the accuracy further in¬ 
creases, varying as ~ The errors achieved 

for both the £ = 0 and £ = 1 channels are two orders 
of magnitude better than the target chemical accuracy 
of Ikcalmol”^ = 0.0016au per electron. The proposed 
pseudopotential is therefore sufficiently accurate for scat¬ 
tering between both opposite- and same-spin electrons in 
this two-body system. 


HEG WITH DIFFUSION MONTE CARLO 

The pseudopotential was calibrated using the exactly 
soluble two-body scattering problem and tested against 
the analytic solution of two electrons in a parabolic trap. 
We now study a system that cannot be solved analyti¬ 
cally: the HEG. We must rely on a numerical approach to 
determine the ground state energy, allowing us to expose 
the computational benefits of using a pseudopotential. 
We first study the HEG with DMC as this is the leading 
approach for accurate calculations of the ground state 
energy [HHUiT] . 

We have used the GASINO quantum Monte Carlo code 
[42] to perform variational and diffusion Monte Carlo 
(VMC and DMC) calculations [U [38]. The Metropolis 
algorithm is used in the VMC method to generate a set 
of electron configurations distributed according to the 
square modulus of the trial wavefunction over which the 
local energy is averaged. In the DMC method, an initial 
wavefunction is evolved in imaginary time, which projects 
out the ground state. The antisymmetry of the wavefunc¬ 
tion is imposed via the variational fixed-node approxi¬ 
mation, in which the nodal surface remains unchanged 
during the evolution. The simulation proceeds with con¬ 
figurations undergoing drift, diffusion and birth/death 
processes, which simulate the evolution of the wavefunc¬ 
tion in imaginary time. DMC provides an upper bound 
on the energy that is lower than the VMC bound calcu¬ 
lated with the same trial state. 

We focus on a three-dimensional homogeneous elec¬ 
tron gas with = 57 electrons and density 

?'s = ?'o/aB = 2 with Ob the electron Bohr radius. The 
calculation is performed in a periodically repeated simu¬ 
lation cell and the interaction energy is calculated us¬ 
ing Ewald summation [431 Slj- We first construct a 
variational wavefunction ip = D that is the prod¬ 
uct of a Jastrow factor e*^ and a Slater determinant 
D = where A is the 

anti-symmetrization operator that accounts for fermion 
statistics. The lowest energy plane-wave states k are used 
to form the orbitals and periodic boundary conditions are 


(a) Time-step error (b) Error with cutoff radius 
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FIG. 3. (Color online) (a) Upper: The error in the energy 
of the HEG with DMC time step. The y-axis origin is set 
by the extrapolation to zero time step of the energy obtained 
with the Coulomb interaction. The red error bars show the 
Coulomb interaction and the blue error bars the pseudopo¬ 
tential, the solid lines show linear extrapolation to zero time 
step. Lower: The uncertainty in energy predictions, with the 
lines showing the fit. (b) The error in the energy with 

cutoff radius for different time steps r G {0,...,!}. Each 
curve has a minimum with cutoff radius, the locus of these 
minima with varying time step is tracked by the green dashed 
line, (c) The relative statistical uncertainty with cutoff ra¬ 
dius in DMC. (d) The spin-resolved pair correlation functions 
for same-spin ((?-|-|-) and opposite-spin (g^t) electrons for the 
Coulomb interaction are shown in red and those for the pseu¬ 
dopotential are in blue. The dotted blue curve shows the an¬ 
alytic correction (g^J’°'^^) applied to the cyan gtlPP obtained 
directly from the pseudopotential, (e) The error in the energy 
and speedup obtained with density. 
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which includes strongly repulsive electron-electron cor¬ 
relations. We describe J by a polynomial expansion of 
order = 8 in the electron-electron separation, [15] 
and Lu is a cutoff length. The behavior of the Jastrow 
factor at electron coalescence can be fixed by the Kato 
cusp conditions m-, for the Coulomb potential we can 
remove the cusp by setting uiap = 3 uoa /3 + 1/2 for an¬ 
tiparallel spins (a ^ /3) and uiaa = Srtoact + 1/4 for 
parallel spins. However, this scheme leaves a remnant 
discontinuity in the local energy. On the other hand, 
the pseudopotential is smooth at r = 0, so there we set 
uia /3 = 8 uoai 3 to ensure that the wavefunction is smooth 
at electron coalescence. The higher order terms in the 
Jastrow factor {Mi> 2 ,a/ 3 } provide the freedom to account 
for longer-ranged correlations. We also add backflow cor¬ 
relations in the Slater determinants using the substitu¬ 
tion r, Ti + - rj|)(ri - r^) with 

V{r) = (1 - r/Lnfe{Lr, - r)J 2 kZo^kr^^ where is a 
cutoff length, and the expansion in variational param¬ 
eters Ck is up to order = 8 [ISj. The variational 
coefficients {wfcct/ 3 , c^, Lu, are optimized using VMC 

[461147]. 

The VMC wavefunction was used as the trial state for 
the DMC calculation. DMC propagates the electrons in 
time step increments r governed by Schrodinger’s equa¬ 
tion in imaginary time. The evolution with a Coulomb 
interaction must have a small time step to properly sam¬ 
ple the rapidly changing local energy near the electron 
cusp [12]. All DMC calculations were performed with at 
least 1000 walkers. We use the percentage of the correla¬ 
tion energy Eq retrieved as the measure of the accuracy. 
There are two main sources of error, firstly the underlying 
VMC trial wavefunction is not exact, having a variance 
in the local energy that introduces a 

systematic error in the DMC estimate of the ground state 
energy of AE = aai^r [iS], where a is a system dependent 
constant. Secondly, because DMC follows a random walk 
there is a statistical uncertainty cte = 
where 5 is a system dependent constant, that can be re¬ 
duced by taking more samples N. Both sources of error 
increase with the variance of the local energy |49j , which 
for the Coulomb potential is ctl = 7.6 x and for 

the pseudopotential (with c = tq) is ctl = 2.4 x 10“^i?c- 
Using the pseudopotential has resulted in a drop in ctl 
by a factor of 3.2, which should reduce both the system¬ 
atic and statistical errors. To expose this we now vary 
another parameter that enters both sources of error: the 
time step. 

In the upper panel of Fig. [^a) we first examine 
the systematic error in the energy. The extrapo¬ 
lates of the ground state energy to zero time step 
for the Coulomb and pseudopotential interactions agree 
to within 0.013%Ac = 0.0012au per electron [15] . 
This is better than our goal of chemical accuracy of 
IkcalmoU^ = 0.0016 au per electron. Calculation with 


the Coulomb interaction and pseudopotential both have 
the expected AE = aanT linear variation of energy with 
time step, though the slope for the Coulomb interaction 
is 3.5-times as steep as for the pseudopotential interac¬ 
tion. This is consistent with the Coulomb interaction 
having a ctl that is 3.2-times as large. Now that we have 
confirmed the analytical form for the systematic error in 
the energy, we examine the statistical uncertainty that is 
expected to be cte = bui ^/The lower panel of 
Fig. da) confirms that the statistical error is well-fitted 
by a power law, and that the ratio of the fitting 

coefficients is 3.3, consistent with the expected ratio from 
the local energy of 3.2. 

With the behavior of both the systematic and statisti¬ 
cal errors verified, we determine the acceleration offered 
by the pseudopotential. Considering only the statistical 
error, cte = &crL/(T^/^iV^/^), to achieve a target final un¬ 
certainty requires a computational effort that scales with 
the number of samples as N ^ The local energy 
calculated with the pseudopotential has an error of ctl, 
which is 3.2-times smaller than for the Coulomb inter¬ 
action, resulting in a 10-times speedup. However, when 
using the pseudopotential the systematic error is also re¬ 
duced, allowing the calculation to be performed at larger 
time steps, which will also reduce the statistical error 
as ~ \Ye consider these effects on an even foot¬ 

ing by combining the systematic and statistical errors in 
quadrature to give a total expected error AEtot in the 
estimate of the energy of 

AEl, = AE^ + 4 

= ( 2 , 


The systematic contribution to the total error grows with 
time step while the statistical uncertainty diverges with 
decreasing time step. The best compromise between the 
two can be found by minimizing the error with respect 
to time step r to yield 


min(AF;tot) 


31 / 262/3 an 
2^/ea^/3 jyi/s ' 


(3) 


If we aim for a particular target total error the com¬ 
putational effort scales with the number of samples as 
TV = (Tl- TFe pseudopotential reduces ctl by a factor of 
3.2, and therefore the pseudopotential offers a ^ 30 fold 
reduction in computational cost while delivering chemical 
accuracy. 

With the benefits of the pseudopotential established, 
in Fig. I^b) we investigate tuning of the pseudopoten¬ 
tial cutoff radius. Starting with a small cutoff radius, 
the energy has a minimal systematic error at small time 
steps, but the calculation with the Coulomb interaction 
suffers from a large local energy variance and the error 
grows rapidly with time step. As the cutoff is increased 
the variance in the local energy is reduced and the fi¬ 
nite time step error falls until it is minimal at c « tq. 
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At large cutoff radii c > ro the interaction potential is 
insufficiently accurate to reproduce the correct ground 
state energy in the zero-time-step limit. There is now a 
high probability that three electrons will be found within 
the cutoff radius, whereas the pseudopotential was cal¬ 
ibrated for two-body physics. The error therefore in¬ 
creases rapidly, independently of the time-step adopted. 
When selecting the cutoff radius one should also consider 
the impact of the variance in the local energy on the sta¬ 
tistical uncertainty in the final result. In Fig. HJc) we 
show that with increasing cutoff radius the increasingly 
smooth pseudopotential leads to a reduction in the rela¬ 
tive uncertainty. At c = tq the relative uncertainty has 
fallen by the same factor of ~ 3.2 as shown in Fig. |^a). 

In Fig. [^d) we study the modification of the pair cor¬ 
relation function arising from the use of the pseudopo¬ 
tential. The same-spin pair correlation function from 
the Coulomb interaction and the pseudopotential agree 
within 0.5%. The opposite-spin correlation functions are 
identical at separations r > c where the underlying po¬ 
tentials are identical. At r c the pseudopotential is 
smaller than the Coulomb potential, and therefore the 
corresponding pair correlation function is larger. How¬ 
ever, at small separations two-body physics dominates, 
and we can separately calculate the pair correlation func¬ 
tion by solving the same two-body scattering problem 
that we used to form the original pseudopotential. This 
two-body solution can be used to correct the many-body 
estimate of the pair correlation function for the incorrect 
two-body effects, bringing it into agreement with the so¬ 
lution for the Coulomb potential to within 1%. Any fur¬ 
ther deviation can be ascribed to three- and higher-body 
physics that occurs for r < c, which is rare as the elec¬ 
trons are simultaneously Pauli blocked and repelled by 
the strong Coulomb repulsion. 

Having confirmed the utility, robustness, and accuracy 
of the pseudopotential for the electron gas with Ts = 2 
we study the accuracy of the pseudopotential for electron 
gases with densities in the range 1 < fs < 16. With 
the cutoff radius at each density set according to c = 
ro = o-Brs, we compare the ground state energy from the 
pseudopotential with that of the Coulomb interaction. In 
Fig-i e) we see that the pseudopotential is able to deliver 
ground state energies to better than chemical accuracy 
with a speedup by a factor of ~ 30 across a broad range 
of densities. 


HEG WITH CONFIGURATION INTERACTION 

The success of the pseudopotential for studying the 
HEG with DMC motivates us to consider a second com¬ 
plementary approach to examine the HEG, Configuration 
Interaction Doubles (CID) [^[3]. We adopt a plane-wave 
basis for our CID calculations, which offers a robust test 
of the portability of the pseudopotential. CID theory 


(a) Coulomb wavefunction (b) Pseudopot. wavefunction 
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(c) Wavefunction convergence (d) Energy convergence 



FIG. 4. (Color online) (a,b) The CID wavefunction for 
opposite-spin electrons passing through coalescence for the 
HEG at Ta = 2 with increasing plane-wave orbital basis sets 
of size M. The exact solution is shown in red (blue) for 
the Coulomb interaction potential and the pseudopotential, 
(c) The average relative error in the wavefunction, and (d) 
the percentage error in the energy with basis set size for the 
Coulomb potential (red) and pseudopotential (blue) with dot¬ 
ted trend lines. 


starts from the Hartree-Fock ground state and includes 
electron correlations through double excitations into the 
unoccupied (plane-wave) orbitals. In the Coulomb po¬ 
tential, the wavefunction has a gradient discontinuity at 
electron-electron coalescence that must be described by 
a large number M of plane-wave basis states with a com¬ 
putational cost that scales as 0{M^). However, the pseu¬ 
dopotential removes the electron-electron cusp rendering 
the wavefunction smooth, which therefore should require 
fewer plane waves to describe the ground state and in 
turn reduce the computational expense. 

The major computational gain offered by the pseu¬ 
dopotential is to aid the description of the behavior at 
electron coalescence, and therefore we first examine how 
the wavefunction at coalescence of two opposite-spin elec¬ 
trons evolves with the size of the plane-wave basis set. In 
the presence of the Coulomb interaction we compare the 
exact relative wavefunction with that from a finite ba¬ 
sis set in Fig.|^a). The Hartree-Fock wavefunction does 
not include opposite-spin correlations, and therefore the 
relative wavefunction is constant at electron coalescence. 
The description of the gradient discontinuity in the wave- 
function at coalescence improves with increasing basis set 
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size. In Fig. Bb) we repeat the exercise in the presence of 
the pseudopotential. The pseudopotential has zero gradi¬ 
ent at r = 0 so the exact wavefunction is now smooth at 
electron coalescence. This allows the shape of the wave- 
function to be described accurately by a relatively small 
basis set. To quantify the change in wavefunction with 
basis set size, we examine in Fig.|^c) the relative error in 
the wavefunction, spatially averaged within the exchange 
correlation hole, k^r < tt, using 


Atp 

t/j 


Jo 



i’ooir)] 


2 

dTrr^dr, 


(4) 


where is the relative wavefunction on coalescence 

of two opposite-spin electrons with separation r, calcu¬ 
lated with CID and a basis set size M. The average 
error for the Coulomb potential falls slowly with increas¬ 
ing basis set size. However, a proper description of the 
wavefunction at coalescence requires a plane-wave basis 
set with a wave vector of at least ^ l.Sfcp, corresponding 
to a basis set size of M = 57. Fig. [^c) shows that here 
the error in the wavefunction drops markedly and the 
wavefunction is over ten times more accurate than that 
for the Coulomb interaction at the same basis set size. 
With large basis sets the wavefunction obtained with the 
pseudopotential converges more rapidly than that for the 
Coulomb interaction. 

Now that we have shown that the pseudopotential fa¬ 
cilitates CID calculations of the wavefunction we study 
the impact on evaluating the ground state energy. Both 
estimates tend towards the same ground state energy, 
confirming the accuracy of the pseudopotential. In 
Fig. I^d) we show that the error in the ground state en¬ 
ergy calculated with the Coulomb interaction scales as 
1/M [in] whereas with the pseudopotential it scales as 
1/M^/^, which is the same improvement as seen with ex- 
plicitely correlated methods m- The pseudopotential 
delivers benchmark chemical accuracy of 0.017%i?c = 
0.0016 au per electron with a ^ 50% smaller basis set and, 
since the computational cost of CID scales as 
this corresponds to a speed-up of a factor of ~ 32. Even 
greater computational gains could be expected at higher 
levels of target accuracy. 

The pseudopotential has contributed to reducing the 
basis set size required in a CID calculation. This benefit 
is expected to be carried over to more accurate configura¬ 
tion interaction approaches, for example coupled cluster 
that overcomes the errors introduced into CID by un¬ 
linked diagrams m- Here we adopted a plane-wave basis 
set, however applications of configuration interaction to 
molecules often express the wavefunction in a coordinate 
basis set centered on the atoms. The pseudopotential 
takes a smooth polynomial form so the two-electron in¬ 
tegrals could be evaluated efficiently as summations over 
the Boys function ISTj. 


(a) Atom energy 
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FIG. 5. (Color online) (a) Error AE in the total energy with 
cutoff radius. The red points are for the isolated Li atom and 
blue are for the Be atom. The gray shading denotes where the 
results attain chemical accuracy, (b) The error in the total 
energy of the Li (red) and Be (blue) ions, (c) The error in the 
ionization energy of a Li (red) and Be (blue) atom, (d) The 
speedup of the DMC calculation. The red points show the Li 
atom and blue the Be atom. 


LITHIUM & BERYLLIUM ATOMS 


The HEG is arguably the most important model of in¬ 
teracting electrons. However, in real systems, the back¬ 
ground charge density due to the atomic nuclei is non- 
uniform and so the electron density varies in space. In 
order to study the performance of the pseudopotential in 
an inhomogeneous system, we perform DMC calculations 
of the energy of the lithium and beryllium atoms. These 
atoms are simple real-life systems that could expose er¬ 
rors introduced by three-body scattering. Accurate ref¬ 
erence results from analytic integration and recursion re¬ 
lations |52l |53| are also available, making these systems 
an ideal test bed for evaluating the performance of the 
electron-electron pseudopotential. 

The trial wavefunction is constructed from single¬ 
particle orbitals in a Gaussian basis set generated by an 
all-electron calculation performed using CRYSTAL [50] . 
The trial wavefunction consists of a determinant of DFT 
orbitals multiplied by a Jastrow correlation factor. The 
parameters in the Jastrow factor are optimized using 
a variance minimization technique |15j . The optimized 
VMC wavefunction is used as a starting point for a DMC 
calculation. 







We first study a solitary Li atom, containing one down- 
spin and two up-spin electrons. We present our esti¬ 
mates for the ground state energy in electronvolts for 
ready comparison with the real-life system. In Fig. [^a) 
we show the variation of the accuracy compared with 
the pure Coulomb interaction. The energy for the ex¬ 
act Coulomb system, —203.379 eV, agrees with refer¬ 
ence results from analytic integration and recursion re¬ 
lations [SH 123] within 0.024 eV per atom. The error de¬ 
creases as the cutoff radius is reduced. If we aim for an 
error of order chemical accuracy (0.025 eV per atom) we 
require c < 0.03ao. Fig. [^d) shows that relative to the 
calculation with the Coulomb interaction, the smoother 
pseudopotential reduces the local variance and therefore 
accelerates the calculation, with greater effect for larger 
cutoff radii. The pseudopotential offers a speedup by a 
factor of ^ 5 while still attaining chemical accuracy. 

The results for the Be atom follow the same trend as for 
the Li atom. For the Be atom we predict a ground state 
energy of —398.932 eV per atom, again within 0.020 eV 
of reference results from analytic integration and recur¬ 
sion relations [SU [53]. The pseudopotential performs 
slightly better for the Be than the Li atom, possibly due 
to the increased prevalence of electron-electron relative 
to electron-ion interaction terms. We also determine the 
energy of the Li+ and Be+ ions in Fig. [^b). The er¬ 
ror is now significantly reduced due to the removal of 
the three-body error for Li’*', and its reduction for Be"*". 
The growth of the error in the energy estimate is sim¬ 
ilar to that for the Li and Be atoms. This means that 
in Fig. j^c) the magnitude of the error in the ionization 
energy grows with cutoff radius. We attain chemical ac¬ 
curacy (0.025 eV per atom) at c < O.OSoq. 

For a fixed target accuracy the speedup of the pseu¬ 
dopotential calculation for the Li and Be atoms is smaller 
than for the HEG. This is because in the HEG we fo¬ 
cused on the error per electron, whereas here we focus 
on the error per atom, which includes three or four elec¬ 
trons, therefore inflating the error. However, even if we 
ignore this, the electron-electron pseudopotential offers 
a 5-times acceleration for high accuracy work, whereas 
for example, for high throughput structure prediction 
calculations an order of magnitude less accuracy is re¬ 
quired [53] so a pseudopotential would offer a 50-times 
speedup. For a molecule chemical accuracy typically re¬ 
lates to the energy difference between two configurations 
rather than total energy for which the pseudopotential is 
expected to be more accurate. 

DISCUSSION 

We have developed a pseudopotential for the repul¬ 
sive Goulomb interaction. The pseudopotential delivers 
accurate scattering states for incident wave vectors and 
angular momentum channels found in an electron gas. 


while its smoothness accelerates computation. With the 
cutoff radius set to the typical electron separation the 
pseudopotential delivers the correct many-body physics, 
and within the cutoff radius two-body physics dominates 
where predictions for the exchange correlation hole can 
be corrected analytically. The cutoff radius can be re¬ 
duced to zero, making the pseudopotential systematically 
improvable. The pseudopotential was shown to deliver 
chemical accuracy for the HEG and to accelerate both the 
DMG and GID methods by a factor of ~ 30. The pseu¬ 
dopotentials were also shown to accelerate the calculation 
of the isolated lithium and beryllium atom by a factor of 
5 for high accuracy work, and in situations where lower 
accuracy is required, for example high throughput struc¬ 
ture prediction calculations, the pseudopotentials would 
provide a 50-times acceleration. 

The performance and simplicity of the electron- 
electron pseudopotential makes it portable across many- 
body techniques such as VMG, DMG, truncated GI, cou¬ 
pled cluster theory, and Mpller-Plesset theory. The for¬ 
malism developed can be applied more widely in scat¬ 
tering problems in condensed matter to develop pseu¬ 
dopotentials for dipolar interactions and also the contact 
interactions found in atomic gases [25]. The approach 
can also be applied to classical physics, for example the 
Coulomb interaction studied here has the same force law 
as Newtonian gravity used in simulations of galactic dy¬ 
namics |55| . Here a pseudopotential could overcome the 
high computational cost and correctly capture the mo¬ 
tion of stars during close encounters. 

The authors thank George Booth and Pascal Bugnion 
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College. 
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